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1 . INTRODUCTION 


In anticipation of long-duration, manned space flights, NASA has 
been investigating Closed Ecological Life Support Systems (CELSS). The 
CELSS would contain a combination of biological, chemical, and mechani- 
cal components. It would provide for the crew's nutritional, atmos- 
pheric, and waste processing needs. The system is closed to mass but 
open to energy. Its goal is the long-term survival of the human crew. 

The CELSS can be thought of as containing a replenishing food sup- 
ply (usually plants), a waste processor, the human crew, and various 
storage tanks. In addition to supplying food for the crew, the plants, 
through photosynthesis, remove carbon dioxide and provide oxygen to the 
atmosphere. As the humans consume the food, they partially reoxidize it, 
using oxygen from the atmosphere and generating carbon dioxide. The 
waste from the humans is run through a waste processor, completing the 
oxidation that the humans started. In this way, the human/waste proces- 
sor components complement the growing food components. The mass flow 
follows a loop and the mass is conserved. 

After considering this general view of a CELSS, the following ques- 
tion was raised: Where might problems develop in system behavior due to 
these component relationships? There are three areas of concern. First, 
how do the long time delays of plant growth affect the ability to con- 
trol the system and insure its survival? Second, what are the effects of 
nonlinearities on the system behavior? Lastly, does system mass, and its 
relation to storage tank capacities, determine the possible system 


behaviors? 
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To address these questions, a series of abstract dynamic models of 
a CELSS were developed [2], They were used to investigate the interac- 
tion of long-term dynamics with finite size storage tanks. The variety 
of missions a CELSS can be used on increases as the total system mass 
and size is reduced [5]. However, the smaller the storage tanks, the 
greater the possibility of a tank overflow during a component failure. 
The series of CELSS models looks at the dynamic consequences of such a 
component failure as a function of tank capacity and control scheme. 

Since no CELSS has been constructed, and the models' nonlinear 
state equations do not resemble any common engineering system, an inves- 
tigation of the abstract models was undertaken. In the simplest model 
there are 5 state variables. In addition to nonlinear functions for 
plant growth, there are 7 switching functions representing empty and 
overflowing storage tanks. While the state space is bounded, the deter- 
mination of the existence of equilibrium points is quite time consuming. 
This simple 5 dimension system only permits one harvest per growing 
time. The more realistic models have multiple harvests in a growing 
period, resulting in 8 to 182 state variables. Even in the case of only 
8, the traditional approach to locating equilibrium points is not prac- 
tical. 


To explore all possible behaviors of these systems, a Monte Carlo 
search of admissible initial conditions was performed. This uncovered 
all possible equilibrium points (and some more complex attracting 
regions) and showed a pattern of bifurcations and reverse bifurcations 
as the storage tank capacities were varied. The relationship between 
system behavior and tank capacity could now be shown. 
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It is important to be able to tell when one controller is superior 
to another from the viewpoint of system survival. Assuming that the 
major threat to a CELSS is a component failure, and that a failure moves 
the system along a random vector in state space, the best controller is 
the one that recovers from the largest variety of failures. This is 
equivalent to returning the system to the desired equilibrium from the 
largest region of the state space. The component failure is viewed as an 
event that moves the system away from the desired behavior. It is 
assumed that the component can be repaired and that it is then the 
controller's job to return the system to its desired operation. 

The domain of attraction of an equilibrium point is the region in 
state space that is attracted to the equilibrium point. Initial condi- 
tions within this region tend toward the equilibrium point asymptoti- 
cally. The attractors do not have to be points, they may be limit 
cycles, or higher-dimension attractors. The present use of domain of 
attraction is of interest only in nonlinear systems because in linear 
systems it is trivial, either being the entire state space or null. 

The domain of attraction can be refined so that a settling time, 
overshoot, or other criteria are included. The resulting region is 
called a domain of performance. While the domain of attraction is not 
particularly enlightening when used with linear systems, the domain of 
performance is useful in investigations of both linear and nonlinear 
systems. 

The control-design goal for the CELSS may be stated as developing 
the controller that gives the system the largest domain of attraction 
around the desired operating point. If this includes the entire 
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admissible state space, then system survival is guaranteed for all pos- 
sible component failures. 

This dissertation explores the characterization and use of the 
domain of attraction and its subset, the domain of performance. Investi- 
gations of the domain of attraction are performed by Monte Carlo 
searches of admissible initial conditions, identifying those that lead 
toward the desired equilibrium point (pass) and those that do not 
(fail). The resulting pass region can then be measured (total volume, 
minimum radius, maximum radius) and a measure of interest can be optim- 
ized. This technique is not greatly burdened by high dimensional systems 
because a random search is used and the performance measure is a set of 
scalars. 

It is possible to mathematically calculate the error associated 
with the Monte Carlo determination of the domain's volume when the sam- 
ple distribution is known. If the domain has a regular shape and low 
dimension, the volume can be used to obtain a minimum and maximum 
radius. The minimum radius is a reasonable measure of the system's abil- 
ity to recover from random perturbations of the state due to a component 
failure. The maximum gives some measure of the elongation of the 
domain. 

For the CELSS models the domain of attraction has an irregular 
shape and a high dimension. The constant mass constraint results in a 
nonuniform sample distribution of initial conditions that is analyti- 
cally unknown. Hence, direct mathematical determinations of volume, 
minimum radius, and associated errors are impossible. Local densities 


are used to obtain volume information from the counts of initial 
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conditions that are contained in the domain. Similarly, the minimum 
radius is found through Monte Carlo simulations and its error is esta- 
blished statistically. 

In the following chapter, the domain of attraction is compared with 
more conventional measures of system stability and performance. Methods 
for measuring the domain and its associated error are shown for both 
uniform and nonuniform distributions. In Chapter 3, controller design 
using the domain as a performance measure is contrasted with more con- 
ventional control designs. Finally, in Chapters 4 and 5, these tech- 
niques are used to examine the CELSS models' behavior as well as some 
properties of CELSS control. 
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2. THE DOMAIN OF ATTRACTION 
i.*JL Stability , Performance , and The Domain 

A system’s equilibrium point (or more complex attractor) is stable 
if nearby points approach it in some specified fashion (e.g.: Lyapunov 
stable, asymptotically stable, etc.). The domain of attraction is the 
entire region about the stable equilibrium point that tends toward it. 
In linear systems, a stable equilibrium point attracts the entire state 
space. Nonlinear systems may have many attractors, each with their own 
domain of attraction. In these cases, both the local and global stabil- 
ity properties are important. 

Consider a system that has its state perturbed suddenly from a 
stable equilibrium. If this is a linear system, it will eventually 
return to the equilibrium. However, if the system is nonlinear, the 
perturbed state may be outside the domain of attraction of the original 
equilibrium. Then the system would not return to the starting equili- 
brium point. Therefore, the size of the domain of attraction in a non- 
linear system is a measure of its ability to recover from displacements 
in its state. 

By including performance criteria, a domain of performance can be 
found. This is a subset of the domain of attraction. The size of the 
domain of performance is an indicator of both linear and nonlinear sys- 
tem performance under impulse perturbations of the state. This is 
because many performance criteria are inherently nonlinear (such as set- 
tling time). Since they do not scale with initial conditions, the 
domain of performance may not be global in stable linear systems. 
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In this dissertation the domains of attraction and performance will 
be used to investigate system behavior. In addition, the measures of 
the domain will be used as a selection parameter in controller design. 

It is very difficult to find necessary and sufficient conditions 
for the stability of nonlinear systems. For example, Lyapunov's second 
method gives only sufficient conditions. Therefore, the region where 
the Lyapunov function exists is a conservative estimate of the domain of 
attraction. Since this is only a sufficient condition, it may be possi- 
ble to find Lyapunov functions that will approximate the domain of 
attraction better. 

Lyapunov's first method can also be used to determine the size of 
the domain for a nonlinear system [17]. The system is decomposed into 
linear and nonlinear parts about the equilibrium point of interest. A 
Lyapunov function is found for the linear part. The system is locally 
asymptotically stable if the gradient along trajectories due to the 
linear and nonlinear parts is negative definite. The domain of attrac- 
tion is approximated by finding the region where the linear contribution 
to the gradient dominates that of the nonlinear part. Since this is 
based on the second method, it is only a sufficient condition and the 
domain estimate is conservative. 

A method for designing controllers and obtaining Lyapunov functions 
for nonlinear systems has been presented by Su, Meyer, and Hunt [16]. 
The nonlinear system is transformed into a canonical linear form. A 
controller can then be designed using standard, linear methods. The 
equivalent controller for the nonlinear system is found by applying the 
inverse transformation to the linear system's control. 
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A Lyapunov function for the nonlinear system can be made by apply- 
ing the inverse transformation to the linear system's Lyapunov function. 
Since this function showed global stability of the linear system, it 
also shows stability over the entire region for which the transformation 
applies. Within this region the nonlinear system can recover from state 
perturbations. 

There are three limitations to this approach. First, the transfor- 
mation must be continuous and smooth ill. A system with switching points 
or other discontinuities cannot be transformed into a linear system. 
Secondly, the transformation and results only apply to a limited region 
of the state space. The region can only contain one attractor of the 
uncontrolled, nonlinear system. As a result, cases where the controller 
greatly increases the domain of attraction cannot be detected. Finally, 
a direct measure of the size of the attracting neighborhood is not 
obtained; only its lower bounds. This makes selection of control stra- 
tegies more difficult. 

The technique proposed in this dissertation avoids many of the 
usual problems associated with determining the domain of attraction. By 
relying on simulation and Monte Carlo selection of initial conditions, 
the procedure is relatively insensitive to system dimension, nonlineari- 
ties, and complexity of form. By finding all attractors, all possible 
system behaviors are discovered. The domain of attraction about a given 
attractor gives a direct measure of the system's ability to recover from 
state perturbations. Controller design can then be accomplished by 
using a Monte Carlo search of control parameters, effectively creating a 
nested Monte Carlo program. Performance measures are easily included. 
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generating a domain of performance that is a subset of the domain of 
attraction. It is also easy to deal with questions of robustness 
against conditions such as parameter uncertainty or disturbance rejec- 
tion. 


2,2 Measures of the Domain 

The domain of attraction is the set of points in state space that 
has trajectories that approach an equilibrium point (or other attractor) 
asymptotically [9]. The domain of attraction is also referred to as the 
basin of attraction [4, 10, 14] and the domain of stability [91. This 
region can be mapped out through simulation using randomly selected ini- 
tial conditions from the admissible state space. Initial conditions that 
are attracted to the given equilibrium are labeled as passes; those that 
are not attracted are called fails. In this section, a method for find- 
ing equilibria and mathematically calculating the volume of the domain 
of attraction will be presented. This applies to both continuous and 
discrete systems. 

Consider the nonlinear, dynamic system: 

x(t) = f(x(t) ,u(t),z) (2.1) 

where x(t) is the state vector, u(t) is an input or forcing vector, and 
z is a vector of parameters. To perform a simulation of this system, 
u(t), z, and an initial condition for x must be specified. While there 
are no restrictions on the form of u(t), it is often taken as a function 
of x(t), giving the system a closed-loop feedback form. The z vector 
represents parameters of the system that can be selected during the 
design process such as controller gains. 
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The equilibrium points are solutions of: 

f(x(t) ,u(t) ,z) = 0 (2.2) 

Finding these solutions can be very time consuming if the dimension of 
the state vector is large or if the function f has many switching 
points. Instead of solving for the equilibrium points directly, a simu- 
lation approach is used with randomly selected initial conditions of x. 
State trajectories are followed until they become (a) stationary within 
specified tolerances, (b) leave a given region, or (c) remain in a 
bounded volume but do not converge to a point. Case (a) is the result of 
a stable equilibrium point. Case (b) occurs when sample initial condi- 
tions are within the domain of attraction of a point outside the region 
being investigated. Higher dimensional attractors (such as limit cycles) 
appear in case (c). 

All attractors (equilibrium points, limit cycles, strange attrac- 
tors, etc.) can now be found for a given region in the state space. How- 
ever, specifying the exact nature of the attractor is not always simple. 
For equilibrium points the trajectory comes to rest within a small 
volume. Limit cycles repeat their motion and can be recognized by a 
repeating pattern of trajectory points. Higher dimensional attractors 
are much more difficult to classify. For example, motion on a torus 
(which has a finite number of frequencies) and motion on a strange 
attractor (which has an infinite number of frequencies) can only be dis- 
tinguished through their spectra or with a calculation of Lyapunov 
exponents [4, 10], In the examples in this and the following chapter, 
only equilibrium points are encountered. More complex attractors occur 
in Chapter 5. 
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Once the attractors have been found, the domain of attraction for 
each may be examined. For simplicity, we will restrict the discussion to 
the domain of one attractor. A random selection of initial conditions of 
x are grouped into passes and fails depending on whether the trajectory 
goes to the attractor of interest. The set of passes indicates the 
region of attraction. This region can be analytically measured. 

To bring this discussion into clearer focus, consider Figure 1. The 
region of interest R (two-dimensional in this case) is the set of 
points: 

R = {(x,y): a£x£b, c£y£d) (2.3) 

The domain of attraction D of point (x ,y ) is: 

o o 

D = t(x,y): x(t)-»x Q , y(t)-»y Q , as t-»oo,and (x,y) € R} (2.4) 


A Monte Carlo integration [6, 12] can be used to determine the area 
of D. If the initial conditions are uniformly distributed on the region 
R, their probability density function is: 


P (x ’ y) = (b-a) (d-c) 


(2.5) 


The probability P that an initial condition is in D (i.e.: a pass) is: 


area D 


p - af gcl U - ft 

area R ' (b-a) (d-c) 


( 2 . 6 ) 


If N random initial conditions are generated, we denote the number of 
passes (in region D) as N^. P can be estimated by: 


P 



(2.7) 
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Using equations 2.6 and 2.7, the area of D can be approximated by: 

- N n 

area D = A Z A = (b-a)(d-e)--jf (2.8) 

Each of the N trials is independent and a binary selection criteria is 
used. Therefore, these are Bernoulli trials with probability P of a 

A 

pass. A is an unbiased estimator of A because: 

- E(N ) 

E(A) = (b-a) (d-c)« ■ ■ ., p ■ = (b-a)(d-c)P = A (2.9) 


The variance of A is: 

varA = [(b-a) (d-c)] 2 varP = A* [(b-a) (d-c) - A] 

N 

a 

where the variance of P is: 

varP = var(N p / N) = — N ~ - 

var p = A* C(b-a) (d-c) - A] 

N • [ (b-a) (d-c) ] 2 

The standard deviation of A is: 

* * 1 /? 
s.d.(A) = [varA] ' 

The precision of the estimate A is proportional to N 


( 2 . 10 ) 


( 2 . 11 ) 


( 2 . 12 ) 


This technique of estimating area (or volume) and error of the 
domain of attraction requires that the initial conditions are uniformly 
distributed. In systems where constraints cannot be substituted in 

explicitly the uniformity of sample points cannot be achieved. An 
implicit solution for this case is demonstrated in the next section. 
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Since the domain indicates the region of initial conditions that 
returns to the attractor, the minimum radius of the domain can be used 
as a worst-case measure for a system whose state experiences a random 
perturbation. This is particularly true when the domain has a complex 
shape such as branches off of a central region. If the perturbation 
moves the system further from the equilibrium point than the minimum 
radius, there is no guarantee that the system is still within the domain 
of attraction. 

2.3, Example : The Damped Pendulum 

A damped pendulum is used to demonstrate some of the techniques for 
finding attractors and measuring their domains. The system has many 
stable equilibrium points, the domain of attraction about one of which 
is examined here. Using a uniform sample of initial conditions, the area 
and error calculations of the last section are compared with a more gen- 
eral technique. Since the minimum radius is a more versatile measure of 
the domain indicating the worst-case system response from a random Per- 
turbation of the state, a general method of finding it and its associ- 
ated error is shown. These techniques are then extended to include the 
case where a nonuniform sample distribution is used. 

The damped pendulum in Figure 2 is represented by the nonlinear 
state equations: 

x = y (2.13) 

y = -sin(x) - b-y 

where x is angular position, y is angular velocity, and b is the damping 
coefficient. Pendulum length and mass have been selected to give unity 
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coefficients. The value of b is set to 0.5 and a time step of 0.05 is 
used with a fourth-order Runge-Kutta integration. The sample region of 
initial conditions is: 

{ (x,y) : -10 < x < 10, -5 < y < 5} (2.14) 

To discover all attractors that can be reached from this region, a 
series of simulations can be performed using initial conditions randomly 
selected from a uniform distribution. The initial conditions are grouped 
according to the equilibrium point that attracts them. A result for 1000 
initial conditions is shown in Figure 3. Points in this region are 
attracted to 7 stable equilibria. Table 1 gives the distribution of ini- 
tial conditions with respect to the attractors, and some representative 
trajectories are shown in Figure 4. 

The origin is selected as the attractor for the investigation of 
the domain of attraction. As a first check of the reliability of the 
Monte Carlo integration, Figure 5 shows that the fraction of initial 
conditions attracted to the origin converges to 31.04%. The total region 
sampled has an area of 200. Therefore, using equation 2.8, the area 
attracted by the origin is 62.08. 

Since the randomly selected initial conditions came from a uniform 
distribution, the accuracy of this area calculation can be established 
using equation 2.12. Area and error results for a variety of sample 
sizes are contained in Table 2. In the standard deviation calculation, 
it is assumed that the true area is 62. 


In most cases, the domain of attraction is not spherical so its 
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Figure 4: Some Pendulum Trajectories 
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volume is not a good measure of system recovery from a random distur- 
bance. The minimum radius is more appropriate because it reflects the 
size of the region where recovery is guaranteed, independent of the 
orientation of the disturbance vector. For systems where the perturba- 
tion has preferred directions, more detailed geometric information on 
the domain may be useful. 

It is possible to get an estimate of the pendulum domain's minimum 
radius by taking advantage of the domain's regular shape. In Figure 6 a 
parallelogram is superposed on the pattern of pass and fail initial con- 
ditions. Its height (y direction) is fixed at 10. The area of the 
parallelogram is: 


A = B-H (2. 15) 

where B is the length of the base (x direction) and H is the height (y 
direction). For a fixed height, errors in area result from errors in the 
base measurement: 


A + a = B-H + b-H (2.16) 

where a and b are the errors in area and base, respectively. 

Using the information on domain area and error in Table 2 with 
equation 2.16, a table of base length and error can be made. The base 
length, however, is not directly a measure of the minimum radius. Using 
the geometry in Figure 7 a measure of the minimum radius r of the domain 
can be found: 

_ _ B*costf 
r 2 


(2.17) 
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Using the error of the area (Table 2), the minimum radius error can be 
calculated: 


a -003(1 
^r = 2H 


(2. 18) 


The minimum radius and error using equations 2.17, 2.18 and Table 2 
for 1000 sample points is: 

r = 1.82 + 0.09 (2. 19) 

Figure 6 was used to set to 53 degrees. 

The above is an approximate method for finding the minimum radius 
relying on the sample distribution being uniform (to find the area) and 
the region being somewhat regular (to relate area to minimum radius). 
There are two problems with this. First, the domain of attraction is 
rarely regular, and second, in higher dimensions a relationship between 
area and minimum radius cannot usually be found. 

To determine the domain's area, minimum radius, and their associ- 
ated errors more directly, the simulation with 1000 uniformly distri- 
buted initial conditions was repeated many times with independent sample 
sets. The area results of these repetitions were stored. The mean and 
standard deviation of the set of area calculations converged before 100 
repetitions were completed (a total of 1000 repetitions were performed). 
The result for 100 repetitions is: 

A 

A = 62.08 + 2.42 (2.20) 

where the number of sample points that passed has been converted into an 
area (equation 2.8). This area agrees with the Table 2 result for 
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100,000 samples (because this is a 100,000 sample result). The error 
agrees with the Table 2 result for 1000 samples. Therefore, the area and 
error calculations using equations 2.8 and 2.12 can be duplicated 
through repeated simulations of the domain of attraction. 

Since the domain’s boundary, and therefore its minimum radius, lie 
somewhere between the pass and fail points, it is not possible to get a 
direct measurement of the minimum radius using Monte Carlo simulations. 
However, a first approximation can be found by locating the initial con- 
dition nearest the attractor that is not in its domain (i.e. a point 
from the fail set). For the pendulum example, this point is shown in 
Figure 8 for 1000 randomly selected initial conditions. The minimum 
radius is the distance from the attractor to the boundary of the domain. 
This point, a fail point, is simply the nearest point that is not in the 
domain. The uniformity of the sample distribution can be used to infer 
the location of the boundary from this point. 

For example, 1000 points uniformly distributed on an area of 200 
gives an expected area per point of 0.2. If this area is in the shape of 
a circle, its radius is 0.25. Assuming the domain's boundary evenly 
divides the region between neighboring passes and fails, the boundary 
would be 0.25 closer than the nearest fail point. 

This approximation of the minimum radius can be refined by replac- 
ing the point's area with an ellipse when the region sampled is not a 
square, but a rectangle. Assuming the ratio of the semimajor to semimi- 
nor axis is 2, the axes measurements for an area of 0.2 would be 0.38 
and 0.18, respectively. Using other shapes for the point's area, such as 
squares, hexagons, etc. can be justified. Such refinements are usually 
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Figure 8: Nearest Fail and Minimum Radius for Domain of Attraction 
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less significant than the errors associated with the minimum radius 
measurement. A more precise refinement of the radius could be obtained 
by concentrating sample points near the nearest fail point. 

Repeated simulations of 1000 uniformly distributed initial condi- 
tions to find the minimum radius (nearest fail point minus the circular 
radius about that point) gives: 

r = 2.26 + 0.06 (2.21) 

(see Figure 8). As before, 100 repetitions were used. While the error 
measurement is similar to that in equation 2.19, this radius is notice- 
ably larger because the domain is not exactly a parallelogram. Rather, 
it has curves on its edges and the minimum radius is at a bulge (see 
Figures 6 and 8 ) . 

In many systems the state variables have constraints that cannot be 
substituted in explicitly. Therefore, a uniform distribution of sample 
points cannot be achieved and an implicit solution for measuring the 
domain is needed. When nonuniform distributions are used, the number of 
pass points no longer directly relates to the volume. The volume can be 
found, however, if the sample distribution is known and can be inverted 
or if the sample point counts are converted to volumes by using the 
local density. The minimum radius of the domain of attraction can be 
found as it was for the uniform distribution. The only difference is 
that the volume per point, used to infer the boundary location from the 
nearest fail point, is now dependent on the local density. 

The pendulum example is repeated with the uniform sample distribu- 
tion replaced by a transformed sample: 
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where U is a sample from a uniform distribution over the range 0 to 1 . V 
is a nonuniform sample over the range 0 to 1 . This is used for both x 
and y coordinates of the initial conditions. A distribution histogram 
for V is shown in Figure 9. 


A 1000 point Monte Carlo simulation is shown in Figure 10 using the 
nonuniform distribution. Clearly, the fraction of initial conditions 
attracted to the origin is not proportional to the domain’s volume. How- 
ever, if the state space is divided into moderately sized, identical 
squares, the initial conditions can be assumed to be uniformly distri- 
buted in each. The number of points in each square gives the local den- 
sity. Therefore, the area of the domain of attraction is: 


A 


K 

1 

i=1 


<N A 


d i 


(2.23) 


where (Np)^ is the number of pass points and d^ is the local density of 
points in square i. The number of squares in the state space is K. The 
area of the domain of attraction using equation 2.23 is: 


A = 64.12 + 2.44 


(2.24) 


The nearest fail point is indicated and the local density is 0.25 
area units per point. Converting this to a circular radius brings the 
domain boundary 0.28 closer than the nearest fail. Repeating this simu- 
lation 100 times and storing minimum radii (nearest fail point minus the 
circular radius about that point) gives the following measure: 


r = 2.34 + 0.06 


(2.25) 




'Figure 9: Nonuniform Distribution Histogram 






Figure 10: Nearest'Fail and Minimum Radius for Domain of Attraction 

Nonuniform Sample Density 
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These results agree with those obtained from the uniform distribution 
(equations 2.20 and 2.21). 

The ability to directly calculate an irregular domain's volume, its 
minimum radius, and their associated error permits this technique to be 
useful for a wide variety of systems. For those where uniform samples 
are not possible (such as constant mass systems), the volume determina- 
tion requires that the sample counts be converted using local densities. 
The minimum radius calculation is virtually unchanged. This means that a 
system's ability to recover from state perturbations can be measured 
regardless of its dimension, the uniformity of the initial condition 
sample set, or the complex shape of the domain of attraction. 
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3. CONTROLLER DESIGN USING DOMAIN OF ATTRACTION 
3..J. Example : The Inverted Pendulum 

An inverted pendulum was chosen to investigate controller design 
techniques (see Figure 11). A variety of controllers are then presented 
to contrast the usual design techniques with those based on the domain 
of attraction. 

The pendulum is constrained to motion in a plane. The controller 
input is the force applied to the cart at the pendulum's base and the 
control goal is to return the cart to the origin with the pendulum 
upright. It is assumed that all four state variables are available to 
the controller. 

The full nonlinear state equations for the. inverted pendulum are: 

x 1 = Xg (3.1) 

2 

. m^L (x^) sinx^ + F - m^cosx^* sinx^ 

X 2 = 2 

m 1 + m 2 + n^cos x^ 


2 

. ^ (m.j + n^'g-sinx^ - m^L(X 2 |) cosx^'sinx^ - F-cosx^ 

x 4 = L ‘ 2 

m^ + m 2 ^ cos 

where x 1 is the cart position, x 2 is the cart velocity, x^ is the pendu- 
lum position (radians), and is the pendulum velocity. Control is 

exerted through force F. The acceleration of gravity is represented by 
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g. The following values are used: 

m 1 = 0.1 (3.2) 

m 2 = 1 

L = 1 

g = 10 

These equations are simulated with a fourth-order Runge-Kutta integra- 
tion with a time step of 0.001. 

For some traditional controller designs, a linear model of the sys- 
tem is needed. Linearizing equation 3.1 about the origin gives: 

= x 2 (3.3) 

. F - 

x - — 

2 " m 2 

*3 = x 4 

. (m^ + m 2 )*g*x^ - F 

X 4 = lTSTJ 

In matrix form this would be: 

x = Ax + bF (3.4) 

where x is the vector of state variables and F is the forcing input. 
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The local stability about the origin can be found by examining the 
poles (eigenvalues) of equation 3.4 (with F =0). The poles are: 

s = 0, 0, +[(m 1+m2 )g / Lm 2 ] 1/2 (3.5) 

Substituting the values of equation 3.2 into 3.5 gives: 

s = 0, 0, 3.32, -3.32 (3.6) 

Hence the origin is unstable (a saddle point) when no control is 
applied. This is equivalent to saying that without control the origin's 
domain of attraction is null. The unstable eigenvalue is due to the 
upright pendulum. The two 0 eigenvalues are due to the inertia of the 
cart. 


A further examination of the linear system shows that it is con- 
trollable using only F as an input. Therefore, the linear system can be 
moved to any point in the state .space in a finite time. Also, the sys- 
tem is observable for output vectors that contain x^. In these control 
examples the entire state vector is known so no state observer is 
needed . 

The control goal is to stabilize the nonlinear (and therefore, the 
linear) system about the origin and to insure that this design can 
recover from the largest class of state perturbations. To show the 
effect of a settling time constraint, a performance criteria is set so 
that the return to the origin (stationary with the cart at zero and the 
pendulum upright) must be within 5 units of time. For investigations of 
control behavior, the domain of attraction will be used with the sample 
initial conditions from the region: 
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{(x-j.Xg.x^Xi}): -10 < x^ < 10, X2=0, (3.7) 

-3.14 < x 3 < 3.14, x 4 =0} 

This region is selected to be two dimensional (x 2 = Xjj = 0) for clarity. 
A three or four dimensional region could be used. 

3.2 Minimum Error Squared Control 

A common technique in modern control design is optimal control 
where a performance measure such as transient time or squared error is 
minimized. Optimal controllers that minimize transient time or squared 
error for linear systems can be expressed in closed form [17]. For most 
nonlinear systems, however, analytic solutions cannot be found. The fol- 
lowing controller will be in a linear state feedback form and it is 
assumed that the entire state vector is either available or observable. 
Since the controller form is specified, this restricted problem is not a 
true optimal control. Rather, it is the minimization of a performance 
measure by performing a parameter search. 

To design a controller for the inverted pendulum that minimizes a 
squared error, a Monte Carlo search can be performed on a range of feed- 
back gains. The controller has the feedback form: 

F = k*x (3.8) 

The state x is a column vector, the gain k is a row vector, and the 
applied force F is a scalar. The error cost function is: 

J - (x^ + w-x^)dt (3.9) 


Thus, error accumulates when the cart is not at the origin or the pendu- 
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lum is not upright. The weighting factor w adjusts the emphasis between 
the two errors. This parameter is part of the controller design and its 
value has a strong effect on the resulting gains. 

An initial condition of: 

X 1 = 5 ’ X 2 = °* X 3 = °* 2 ’ x 4 = 0 (3.10) 

is used to find the control gains that minimize the cost function (equa- 
tion 3.9). To insure that the pendulum is quickly moved toward the 
upright position, a weighting factor of w = 1000 is used. Five hundred 
sets of gains were randomly selected from uniform distributions with the 
following ranges: 

f° < k T < 10 °* 0 < k 2 < 100, (3.11) 

0 < k 3 < 400, 0 < k^ < 200} 

The minimum cost function (J = 78.1) was obtained with the gain: 

k = [15.2, 49.2, 305., 64.4] (3.12) 

This gain was used in a Monte Carlo simulation of 2000 initial condi- 
tions randomly sampled from the range in equation 3.7. The domain 

attracted to the origin in t < 5 is shown in Figure 12. The minimum 

radius is: 


?perf = °-35 i °.°6 (3.13) 

This region is not the total domain of attraction. There is a 
larger set of initial conditions that have trajectories that lead to the 
origin. Figure 13 shows the stability region for this control. Here, all 
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Figure 13: Domains of Performance and Attraction for Minimum Error Squared Control 
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initial conditions get to the origin in t < 20. The minimum radius is: 

? stab = 0,81 - °' 03 (3 ‘ 14) 

Figure 14 shows all 2000 initial conditions. 

To distinguish between the various regions in the state space, the 
following terms are introduced. 

Stable set: 

Initial conditions that are attracted to the equilibrium point. 
Unstable set: 

Initial conditions that are not attracted to the equilibrium point 
(complement of the stable set). 

Pass set: 

Initial conditions that satisfy the performance criteria (a subset 
of the stable set). 

Fail set: 

Initial conditions that do not satisfy the performance criteria 
(complement of the pass set). 

Domain of performance: 

Region delineated by the pass set. 

Domain of attraction: 

Region delineated by the stable set. 

Figure 15 shows a typical trajectory from an initial condition in 
the stable region. The cart motion does not meet the settling time cri- 
teria. Trajectories in the pass region are similar but have faster 
responses. A typical unstable trajectory is shown in Figure 16. The 
failure is dramatic. The pendulum starts spinning rapidly and the cart 
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T I me 

Figure 16: An Unstable Trajectory for Minimum Error Squared Control 
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accelerates away from the origin. At this point, the controller is com- 
pletely confused by the nonlinearities of the system. 

Although this control design minimizes the squared error for a 
specific initial condition, the optimization is very local. This gain 
would not be optimal for other initial conditions since the system is 
nonlinear. The design therefore depends on two parameters: the weighting 
factor and the initial condition. There is no a priori way of determin- 
ing an appropriate choice for them or their effect on the system’s glo- 
bal behavior. In this case, the optimization of a local region occurred 
at the expense of system response from a more global set of initial con- 
ditions. 

3,. 3. Minimum Time Control 

An alternative cost function for the optimal control design is 
minimum time. The previous example showed a very small domain of perfor- 
mance but a reasonable size stability region. Hopefully, using time as 
the cost function will enlarge the region that passes the time limited 
performance criteria. Since the controller form is restricted to linear, 
state feedback, this design is not an optimal control. Rather, it is the 
minimization of the system settling time given this controller form. 

Using a state feedback controller (equation 3.8), the cost function 
is now: 


J = t: x. _< 0.01; i = 1 ,2,3,4 


(3. 15) 


where t is the time it takes the system to settle to within the toler- 


ance band about the origin 
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The initial condition in equation 3.10 was used to find the gains 
that minimize the cost function J. Five hundred sets of gains were ran- 
domly selected from uniform distributions over the range in equation 
3.11. The minimum cost function (J = 2.86) was obtained with the gain: 

k = [48.3, 59.5, 352., 91.0] (3.16) 


This gain was used in a Monte Carlo simulation of 2000 initial con- 
ditions randomly sampled from the range in equation 3.7. The domain 
attracted to the origin in t < 5 is shown in Figure 17. The pass region 
is identical to the stability region. In fact, all initial conditions 
surpass the performance criteria and arrive at the origin in t < 3.5. 
The minimum radius is: 


perf 


= r t . = 0.89 + 0.03 
stab — 


(3.17) 


While this control design insured that the selected initial condi- 
tion met the performance criteria, it did not generate a domain of per- 
formance reaching very far beyond it. The sharp edge of the domain 
(t 3.5 for attracted points) is a result of the local optimization. 
Other initial conditions would generate different results. While this 
design meets the performance criteria for a larger set of state pertur- 
bations than the minimum squared error design, the consequences of local 
design decisions (the initial conditions) on the global result (domain 
size) are not easily predicted. 

3.4 Pole Placement Control 


Pole placement control design permits the manipulation of system 
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Figure 17: 2000 Initial Conditions for Minimum Time Control 
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performance through the selection of closed loop poles. The system model 
is linearized so the design is only insured for behavior near the 
equilibrium point. The controller is in a state feedback form and it is 
assumed that the entire state vector is either available or observable. 
For controllable systems, the system/controller poles can be placed 
wherever desired. 

Using the linearized system described in equations 3.3 and 3.4, the 
control is expressed in a state feedback form (equation 3.8). Combining 
equations 3.4 and 3.8 gives: 

• 

x = (A + b-k) x (3.18) 

where A is the system matrix and b is a column vector. The poles (eigen- 
values) of equation 3.18 can be placed arbitrarily because the system 
(A,b) is controllable [3]. Transforming the system to a control canoni- 
cal form aids in finding the gain k that yields the desired poles. 

It is difficult to determine which pole locations will satisfy the 
system design criteria. The poles are selected in a linear environment 
but the system is highly nonlinear. Therefore, the task is to find poles 
that will permit the nonlinear system to recover from the largest range 
of perturbations with t £ 5. 

For instance, if the poles are placed at s = -1, -1, -5, -5 the 

gain vector is: 

k = [2.5, 6, 59.5, 18] (3.19) 

Equations 3.8 and 3.19 show that the slow poles (-1) cause a moderate 
reaction to cart position and velocity errors while the fast poles (-5) 
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cause a stronger reaction to pendulum position and velocity errors. 

Two thousand uniformly distributed initial conditions from the 
range in equation 3.7 were simulated using the gain in equation 3.19. 
The domain of performance (attracted to the origin in t < 5) is shown in 
Figure 18. It has a central region and two separated wings. The wings 
are isolated from the main body of the attractor by a region that is 
stable but does not meet the settling time limit. The total set of ini- 
tial conditions is shown in Figure 19. 

Given the performance criteria of t £ 5, the controller's behavior 
can be measured using the pass region. The region’s wings and long 
extensions along the x axis do not contribute to the system's ability to 
recover from random state perturbations. Rather, perturbation recovery 
is measured with the minimum radius of the domain of attraction: 

Pperf = °- 45 ± °‘°3 (3.20) 
For comparison, the minimum radius of the stability region is: 

p ,tab * ’-20 i°-°3 ‘3-2D 

To increase the system's ability to recover from perturbations, the 
main body of the domain of performance needs to be expanded toward the 
wings, filling in the stable region. Since this region was not satisfy- 
ing the time constraint, it would seem reasonable to use faster poles to 
improve the system behavior. However, increasing the speed of all the 
poles has been found to shrink the region. For example, using poles at s 
= -5, -5, -5, -5 gives a performance minimum radius of 0.31 +0.05. 
Through trial and error, it was found that good results were obtained 




Figure 18: Doma 
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when an appropriate mixture of fast and slow poles were used. There is 
no direct method to determine this mixture. 

It is generally true in linear systems that faster poles give the 
system a larger domain since they imply that the system recovers quickly 
from perturbations. As shown above, this is not necessarily true of non- 
linear systems. The example with poles at s = -1, -1, -5, -5, recovers 
from a larger group of perturbations in t < 5, yet its linear design has 
some slower poles than the case with all poles at -5. This is because 
the linear approximation only holds for behavior near the equilibrium 
point. The domain investigation uses initial conditions far from the 
origin; here the nonlinear properties of the system dominate its 
behavior. Thus, the decision to use faster or slower poles cannot always 
be intuitively based on the linearized behavior. 

3.5 Monte Carlo Search For Linear State Feedback Gains 

The connection, or map, between the placement of poles in a linear- 
ized design and the resulting behavior in the nonlinear system is nonin- 
tuitive. The limitation in the use of linear state feedback control on 
nonlinear systems is in the absence of this map from design parameters 
to system performance. 

A more direct design approach is used here. The goal remains to 
find the linear feedback controller gains that enable the system to 
recover from the largest set of perturbations. However, instead of look- 
ing for linearized system poles or optimal gains for a specific initial 
condition, a direct search of controller gains is made using the domain 
of performance as the selection filter. This global design method 
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searches for systems that can recover from the largest class of state 
perturbations. 

The controller has the feedback form of equation 3.8. A range for 

each element of the k vector is selected and the elements are set with a 

random search using a uniform distribution over the respective ranges. 
This gain vector is used in a domain simulation of the nonlinear system 
with initial conditions selected randomly from the region of state 
space. At the end of each domain simulation, the minimum radius of the 
domain of performance is stored and another gain vector is picked. The 
largest minimum radius found from this nested Monte Carlo procedure 
corresponds to the system that can recover from the largest set of ran- 
domly oriented state perturbations. 

The gains are randomly selected from the set: 

{0 < k < 50, 0 < k 2 < 50, (3.22) 

0 < k 3 < 200, 0 < k^ < 100} 

where each subscripted k is an element of the gain vector. Five hundred 
gain vectors were selected and each gain was simulated using 1000 ini- 
tial conditions from the range in equation 3.7. 

The best gain was found to be: 


k = [7.2, 13.2, 103.8, 39.4] 


(3.23) 


Figure 20 shows the domain of performance for this gain using 2000 ini- 
tial conditions. The minimum radius is: 


perf 


stab 


1.20 + 0.03 


(3.24) 



Figure 20: 2000 Initial Conditions for Linear Feedback Control 
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The minimum radius of the pass and stability regions are identical. 
They differ only where the cart is far from the origin and the pendulum 
is leaning in the same direction as the cart displacement. In these 
cases, the time to get the pendulum balanced and move the cart back to 
the origin is slightly longer than the performance limit. 

Using the gain in equation 3.23, the poles of the closed loop sys- 
tem (equation 3.18) are at s s -22.66, -1.02, -1.26+1. 22 j, where j is 
the unit imaginary number. Hence, the best mixture of poles is for one 
to be very fast and the other three to be rather slow. 

While this design has a larger minimum radius than that found by 
trial and error placement of poles or local optimizations, it may not be 
the controller with the largest possible domain. The range of gains was 
sampled only 500 times. Further sampling could turn up some that would 
perform better. 

3,. 6^ A Nonlinear Controller 

Of the controllers presented in this chapter, the one that 
recovered from the largest class of state perturbations was the linear 
state feedback using a search over a range of gains. The design did not 
depend on parameter selection and the domain of attraction was used to 
determine global behavior. A critical limitation, however, exists in the 
system' s recovery from perturbations in the direction of pendulum angle 
(x^). To try to increase the angular range, a nonlinear controller with 
nonlinear gain for the angle component of the state feedback will be 
used. 


The nonlinear controller for the inverted pendulum uses the 
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following feedback form: 

F = [k 1 x 1 + k 2 x 2 + fCx^x^ + (3.25) 

a^x^: -b £ _< b 

a 2 X 2 + c: x^ >. b 

a 2 x^ - c: x^ _< -b 

The nonlinear gain function f is shown in Figure 21. 

To reduce the dimension of the gain search, the gains for x^ and x 
are fixed at the values found for the best linear gain in equation 3.23. 
The angular gains are from the range: 

{0 < a 1 < 300, 0 < a 2 < 400, 0 < k 4 < 80} (3.26) 

The break point in the nonlinear gain f is set to b = 0.9. The inter- 
cept c is calculated for each selection of a^ and a 2 so that the func- 
tion is continuous. 

Five hundred gains were randomly selected from a uniform distribu- 
tion over the range in equation 3.26. Each was simulated with 1000 ini- 
tial conditions distributed over the range in equation 3.7. The system 
with the largest domain had gains: 

k 1 = 7.2 (3.27) 

k 2 = 13.2 


f(x 3 ) = 


a 1 = 118.7 
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a 2 = 253. 


\ = * 3.5 

The domain of attraction using 2000 initial conditions is shown in Fig- 
ure 22 and the minimum radius is: 


r perf ” r stab 


1.32 + 0.03 


(3.28) 


For this controller, the domain of performance and the domain of attrac- 
tion are identical over the sampled range of initial conditions. 


While this controller is an improvement over the linear one, it is 
probably not optimal. First, the gains on cart position and velocity 
could be included in the Monte Carlo search. Secondly, all the gains 
could be made nonlinear with adjustable break points. Lastly, a state 
space discrete controller, such as those proposed by Young [18], could 
be used instead. This control design example was chosen to show that 
nonlinear controllers can be designed with the same techniques as linear 
ones. 


3.. 7 Summary 

In designing controllers for nonlinear systems, the nonlinearity 
cannot be ignored. Controllers that are designed with a linear approxi- 
mation may or may not be adequate to control the nonlinear system. The 
linearized design only guarantees behavior near the equilibrium. While 
this is sufficient to determine local stability, it is often inadequate 
to evaluate the system's global behavior. Further, controllers that 
optimize a local behavior often do so at the expense of global perfor- 



0 0 <*> 



Figure 22: 2000 Initial Conditions for Nonlinear Control 
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mance. By using the domain of attraction as a selection filter in con- 
troller design, the global behavior of the nonlinear system is directly 
known. The nested Monte Carlo method permits design of linear and non- 
linear controllers that have large domains of attraction and perfor- 
mance. The elimination of design parameters, such as weighting factors 
or initial conditions, further simplifies the design problem. 

A distinction has been made between domain of performance and 
domain of attraction. For systems that are initially unstable, domain of 
attraction may be sufficient as a performance measure. However, stable 
systems, particularly those that are globally attracting, can only be 
analyzed with the domain of performance. The only difference is in the 
goals; the domain of attraction is used to delineate stability and the 
domain of performance indicates the region that meets the performance 
criteria. 

A tabulation of minimum radii for the various controllers presented 
in this chapter is found in Table 3. 


Table 3 


Inverted Pendulum Controller Summary 


Controller 

A 

r 

perf 

A 

r 

stab 

Minimum Error 

0.35 + 0.05 

0.81 + 0.03 

Minimum Time 

0.89 + 0.03 

0.89 + 0.03 

Pole Placement 

0.45 + 0.03 

1.20+0.03 

Linear Search 

1.20 + 0.03 

1.20+0.03 

Nonlinear Search 

1.32 + 0.03 

1.32+0.03 
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4 THE CELSS MODELS 


4_.J_ Overview 

A Closed Ecological Life Support System (CELSS) is usually con- 
ceived as having mass closure but an external supply of energy. It is 
possible to model such a system using conservation of mass equations 
that describe the storage tank behavior. Flows between tanks can be 
treated as controllable variables. Averner [1] used this approach 
balancing the elemental masses (hydrogen, oxygen, nitrogen, etc.) in the 
system. Stahr, et al. [15] developed a model where bulk masses (water, 
carbon dioxide, edible food, etc.) are followed through the system. The 
latter approach will be used here because it lends itself to examina- 
tions of the interplay between storage tank size, system mass, and con- 
trol design. 

Initial design studies for closed life-support systems concentrate 
on the equilibrium requirements for supporting the crew [11], These stu- 
dies give some indication of mass and volume requirements by specifying 
the flows that will be necessary through various processors, and thus 
give some indication of the minimum unit size. However, the life support 
system must be capable of maintaining vital functions during temporary 
failures of some of its components. Extra storage must be provided, pro- 
cessors must have the capability of operating above (or below) their 
equilibrium flows, and total amounts of flowing masses in the system 
must be specified. This part of the design can only be done by consider- 
ing the system's dynamic behavior as none of these parameters enter into 
the static equilibrium calculation. 
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To examine the dynamic interaction of internal system mass and 
storage tank sizes, the model of a CELSS shown in Figure 23 is used. 
This abstract model is a simplification of the true complexity of the 
system. It does, however, contain some essential components of a CELSS: 
constant mass, finite storage tank sizes, and limited processor capaci- 
ties. 


To understand the model better, its behavior at steady state is 
examined. The crew consumes food at a rate of one unit/day (steady 
state) from food storage. With the crop growing to maturity in 60 days, 
the 6 plant chambers produce a harvest every 10 days. This harvest must 
have an edible mass of 10 units for the system to remain at steady 
state. 


Under normal conditions, one third of the harvest is edible. The 
inedible part is placed in waste storage. As food is consumed, the 
crew’s waste also goes into this tank. The waste is reoxidized in the 
waste processor and the resulting nutrients, water, etc., are placed in 
the plant chambers. To insure adequate growth for the steady state, the 
processor flow must be 3 units/day. 

Both food and waste storages have capacities. If a capacity is 
exceeded, the tank's output flow must be increased. The waste processor 
also has a capacity. If its capacity is exceeded, some of the flow is 
bypassed. Clearly overflow conditions can occur given finite storage 
tanks and the possibility of a component failure. 

The flows in this model are mixtures of solids, gases and liquids. 
Thus, the "nutrient" flow refers to nutrients, water, carbon dioxide. 




Figure 23: A CELSS Model 
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and other material needed for plant growth. The "harvest" contains 
excess water, oxygen, edible and inedible plant matter, and other Bypro- 
ducts of plant growth. The proper mixing of the elements is assumed. In 
this model the plant growth depends only on the rate of the nutrient 
flow. 


The plants in each chamber are at different stages of growth to 
allow harvests at 10-day intervals. As a simplification, an independent 
supply of seeds is assumed. The steady state plant growth is shown in 
Figure 24 (top). This curve follows the general behavior of plant 
growth [133. The plant mass reaches 10% of the harvest mass in the first 
20 days from a nutrient flow of 0.1 units/day. The plants grow faster in 
the second 20 days reaching 50K of its total mass. So far no edible mass 
has grown. In the last 20 days, 2/3 of the nutrient flow goes into edi- 
ble growth. Hence, the result is a harvest that is 1/3 edible. 

If the nutrient flow into the plant chamber is not at the steady 
state two things happen [8]. First, the total plant mass does not 
increase as indicated in Figure 24 (top). This effect is due to conser- 
vation of mass. Second, if this occurs during the last 20 days of the 
growth cycle (when the edible part grows) , the percent of the nutrient 
flow that becomes edible is affected (Figure 24, center). The nutrient 
flow into the 6 chambers is always divided to do the least damage to the 
growing plants. 

If the waste processor’s capacity is exceeded, the bypass flow goes 
into the plant chambers (see Figure 23). This overflow accumulates as 
"inert matter" in the plant chambers and does not contribute to the 
plants' growth. If there is inert matter in the plant chamber during the 
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last 20 days (when the edible part of the plant is growing) its growth 
is inhibited (see Figure 24, bottom). The inert matter is removed from 
the plant chamber during harvest, and is sent to the waste storage. 

Control is exerted on the system by adjusting the nutrient flow to 
the plant chambers. Flow regulation decisions may be based on a variety 
of information about the system. For example, the growth of plants 
could be set to compensate for a poor harvest, or the flow of nutrients 
could be adjusted to return the waste storage tank to its steady value. 
Evaluating the effectiveness of these feedback schemes is critical in 
the controller design. This, along with an investigation of the effects 
of storage tank size will be covered in the following sections. 

The State Equations 

Three versions of the CELSS model were considered: one with three 
plant growth chambers, one with six, and the other with twelve. While 
the steady behavior of these three are similar, their differences in 
harvest rates and storage tank needs are considerable and have far 
reaching consequences. 

The three-chamber model looks similar to the one in Figure 23 
except for the number of growth chambers. The plant growth curves are 
the same as those in Figure 24. It has harvests every 20 days whereas 
the six-chamber version is harvested at 10-day intervals. The plant 
growth time is 60 days for both. Similarly, the twelve-chamber model has 
the same form but is harvested every 5 days. 

The CELSS models consist of state equations that represent mass 
conservation for storage tanks and the plant growth chambers. The state 
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equations are of the form: 

x(t) = f(x(t) ,u(t) ,g(t) ) (4.1) 

where x is the vector of masses, u is the control, and g is the vector 
used to represent switching functions associated with harvesting and 
tank capacities. 

For example, the food storage state is: 

x food * ’in - q out * 8 1 <x food- l) (, - 2) 

The flow into the tank (q. ) is the edible part of the harvest. It is a 

in 

delta function because the harvesting is periodic. The output flow 
^ q out^ is the re 3 uired rate of food g° in 8 to the crew (a continuous 
function). The function 8i( x food ) is a switching function that is used 
to represent an empty or overflowing food tank. The waste storage tank 
has a similar state equation except that its switching function g also 
contains information on the waste processor capacity. 

Each plant chamber has three state equations associated with it. 
The first is the mass of the edible part of the plant. It is set to zero 
for any chamber associated with the first 40 days of the plant's growth. 
During the last 20 days of plant growth the state equation is: 

• 

x edible = q in* f 1* f 2 + g 2 (x edible' (4,3) 

The input flow (q^) is the nutrient flow from the waste processor. 
Functions f^ and f^ are the edible growth functions represented in Fig- 
ure 24 (center and bottom, respectively). f 2 is the inert mass state 
variable associated with each chamber. The switching function g 2 
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represents periodic harvesting. 


The second state variable for each chamber is the inedible part of 
the plant. For plants over 40 days old, it grows by using the nutrient 
flow (q in ) that is not being used for edible growth: 


X inedible = q in* C1 “ f 1* f 2 ] + g 3 (x inedible’ t) 


(4.4) 


The functions f^ and f^ are the same as for the edible state variable. 
The switching function g^ represents periodic harvesting. For chambers 
holding plants that are less than 40 days old, the inedible fraction 
uses all of the nutrient flow: 


x inedible = q in 


(4.5) 


The third state variable associated with each chamber represents 
the inert mass. This is the accumulation of unprocessed waste in the 
plant chambers due to the bypass flow. Its state equation is: 

• 

"inert * "in * ( "- 6) 

where q in is the bypass flow sent to the chamber. Since the inert mass 
affects edible growth as an accumulated quantity (see Figure 24, bot- 
tom) , it is most advantageous to divide the bypass flow evenly between 
all plant chambers. The function represents periodic harvesting. 

The total number of state variables depends on the number of plant 
chambers. When three are used, there are 9 mass storage variables (1 
edible, 3 inedible, 3 inert, 1 food storage, 1 waste storage). For the 
six-chamber model there are 16 (2 edible, 6 inedible, 6 inert, 1 food 
storage, 1 waste storage). The twelve-chamber model has 30 (4 edible. 
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12 inedible, 12 inert, 1 food storage, 1 waste storage). Since the sys- 
tem has constant mass, only n-1 of the n mass storage variables can be 
independently selected. These are the state variables. The explicit 
time dependence of the harvest in the continuous-time system appears as 
an added state variable, creating an "augmented state space" [173. This 
augmented state vector has dimension n. 

The harvesting time interval depends on the number of plant 
chambers. With three the interval is 20 days; with six it is 10 days; 
and with twelve it is 5 days. The harvest looks like a periodic delta 
function. 

By sampling the continuous-time system state at the harvesting fre- 
quency, a return map (Poincare surface of section) can be made [10]. 
This map is a discrete-time system where the time dependencies of the 
forcing functions are now implicitly included in the formulation. The 
system (equation 4.1) now appears as: 

x(k+1 ) = x(k) + T*[f d (x(k),u(k),g(k))3 (4.7) 
where k is the time index representing steps of T units. 

The continuous-time system function f is now replaced with an 
equivalent discrete-time function f^. The key difference in the two is 
in how delta functions are treated. In the former case they are 
evaluated at the appropriate time. However, with the larger time steps 
of the discrete-time simulation, the delta functions need to be approxi- 
mated. This is accomplished by evaluating them first in the sequence of 
calculations for each time step. In this way the discrete-time simula- 
tion generates the same behavior as the continuous one in most 
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situations. 

One limitation with the discrete-time calculation is that it cannot 
properly simulate continuous phenomena as observed in the behavior of 
the system during a component failure or tank overflows due to continu- 
ous (not delta function) inputs. In the examples that follow, there are 
no occurrences of such overflows. But when a component failure occurs, 
the continuous-time simulation must be used to follow the system trajec- 
tory. The discrete-time simulation can then be used to examine the sub- 
sequent evolution of the system once the component is repaired. Thus, 
the system's ability to recover from random component failures can be 
investigated through simulations of random initial conditions using the 
discrete-time calculation. 

£.3, Equilibrium Behavior 

The three CELSS models described in the previous sections were 
designed with a steady state that provides harvests equal to the crew 
demand of one unit/day. The relationship between the continuous-time and 
discrete-time simulations of these models can be seen by examining the 
system behavior at this equilibrium. The three- and six-chamber models 
are demonstrated in this section. Similar behavior is observed in the 
twelve-chamber model. 

Figure 25 shows the three-chamber model operating at steady state. 
The continuous simulation's behavior is represented by the plotted line 
and the circles are the output of the discrete simulation. Figure 26 
shows the equivalent steady state for the six-chamber system. 


Each harvest causes a vertical jump in the storage curves while 
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Figure 25: Continuous- and Discrete-Time Equilibrium; 

Three-Chamber Model 
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Figure 26: Continuous- and Discrete-Time Equilibrium; 

Six-Chamber Model 



74 


continuous output flows cause smooth downward slopes. This sawtooth pat- 
tern indicates that the continuous simulation's equilibrium state is not 
a point but a cycle. The discrete simulation has a time step equal to 
the harvest interval so its output does not show this cyclical behavior. 
Its steady state is an equilibrium point. The stability of this steady 
state will be examined in the next chapter. 

In these examples, the food and waste tanks are never empty. This 
buffer level is large enough to maintain the output flows for 10 days. 
This buffer level can be set wherever desired. However, the survivabil- 
ity of the system to survive failures depends a great deal on the size 
and location of these buffers. 

The system mass and size of the storage tanks needed depend on the 
number of plant chambers the model has. Both the three- and six-chamber 
models have food and waste buffers of 10 and 30, respectively. However, 
the mass tied up in the plant chambers is not the same. The three- 
chamber model has 96 mass units in its plant chambers at steady state 
while the six-chamber model has 81 units. Further, the three-chamber 
model's larger harvest mass (20 days of food and inedible waste) 
requires larger food and waste tanks than the six-chamber model does. 
With these static considerations, the lighter and smaller six-chamber 
system would appear to be a better choice. However, in the next chapter 
the dynamic behavior of the two models is compared and the superiority 
of the six-chamber system is not as obvious. 

4.4 A Processor Failure and the Need For Control 


Some of the system dynamics can be seen when a waste processor 
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failure is considered. In this example, the supply of water and 
nutrients to the plants is stopped from the fifth to the fifteenth day. 
During this 10-day failure the output of the waste storage is set to 
zero to avoid a bypass of the processor. From a static design viewpoint 
the steady state contains enough food and waste in their respective 
buffers to ride out the 10-day processor failure. Control is exerted on 
the system through adjustments of the output flow of the waste storage. 
All storage tanks are large enough to prevent overflows. 

One strategy is to maintain the waste storage output flow at its 
steady state value. Using the continuous-time simulation, the effects of 
this no-control strategy on a 10-day failure are shown in Figure 27 for 
the three-chamber model and in Figure 28 for the six-chamber model. The 
system returns to an equilibrium in 60 days. This is not the original 
steady state, however. The food buffer of 10 units is gone and the waste 
buffer has increased from 30 units to 40. This transfer of mass leaves 
the system in a configuration that would be disastrous if another pro- 
cessor failure occurred. 

In addition to moving to a new steady state, the three-chamber 
model needs a waste storage capacity of at least 93 units to absorb the 
transient without causing an overflow. For the six-chamber model a waste 
capacity of 74.25 is needed. This dynamic behavior could not be 
predicted from the static design. 

To return the system to its original steady state with 10 day 
buffers, some sort of control action is needed. Figure 29 shows a feed- 
back control that adjusts the nutrient flow to the plant chamber. If the 
harvest is too small, the controller increases the nutrient flow which 
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Figure 28: Processor Failure with No-Control; 

Six-Chamber Model (Continuous-Time) 
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Figure 29: CELSS with Proportional Control 
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in turn increases plant growth. The gain adjusts the magnitude of the 
change in nutrient flow relative to the harvest error. The feedforward 
of the set point avoids offset errors in the control. 

The feedback controller is expressed mathematically as: 


u 


u steady 


M ♦ k(1 - y/y set H 


(». 8 ) 


where u is the nutrient flow to the plants and y is the edible harvest. 
The proportional gain is k, u stea( jy * s the steady state nutrient flow, 
and y . is the desired edible harvest mass. At each harvest the con- 
troller calculates the nutrient flow which is then maintained until the 
next harvest. Thus, this is a discrete-time controller. If the control 
function results in a nutrient flow large enough to harm the plants, the 
controller applies the flow that is at the maximum point of the growth 
curve (see Figure 24, center). 


This controller is applied to the processor fail example. With a 
gain k = 1 the system behavior can be seen in Figures 30 and 31, for the 
three- and six-chamber models, respectively. The food and waste buffers 
are recovered after the failure, however, the recovery time is between 
100 and 120 days. Such a long delay after a failure is due to the system 
dynamics. Even though the longest time constant in the model is the 
plants' 60-day growth time, the system shows evidence of the disturbance 
for much longer periods. As was the case with the no-control example 
(Figures 27 and 28), larger storage tanks are needed to avoid overflows 
than is indicated by a static design. 
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Figure 30: Processor Failure with Proportional Control; 

Three-Chamber Model (Continuous-Time) 
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Figure 31: Processor Failure with Proportional Control; 

Six-Chamber Model (Continuous-Time) 
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5 DOMAIN OF ATTRACTION IN CELSS ANALYSIS 

The methods of system analysis using the domain of attraction 
presented in Chapters 2 and 3 are particularly suited to investigating 
many aspects of the CELSS models. Their high dimensionality and non- 
linear state equations make them difficult to analyze by any other tech- 
nique. The goal is to investigate the dynamic consequences of a change 
in the waste storage capacity and the total system mass, as well as how 
information is used by the controller. 

A random component failure is simulated as a random initial condi- 
tion. The size and shape of the region of initial conditions that are 
attracted to the specified operating point indicate the system's ability 
to recover from random failures. 

Three-, six-, and twelve-chamber models are considered. The waste 
processor capacity is 4.5 units/day. Waste flows exceeding this amount 
are bypassed. 

5_. 1_ Uniformity of Sample Points 

Since the system has constant mass, the state space is bounded. 
Therefore, not only can the entire state space be sampled in the selec- 
tion of initial conditions, it is also possible to find all attractors 
in the system. The constant mass constraint, however, makes it impossi- 
ble to get a uniform distribution of initial conditions. 

Consider a system with n mass storage variables: x 1 through x n> 
Since the system has constant mass, only n-1 of these are needed to 
describe the state. Therefore, it is possible to randomly select the 
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first n-1 mass variables, but this fixes the value of x . Given the 

n 

system's mass, it is possible to set bounds on the n state variables: 

0 i x i i M i* 1 = 1.2, ... n (5.1) 

where NL is the maximum value for x^. The M i may be simply set to the 
total mass or more refined calculations can be done to set each vari- 
ables' maximum. 

The set of initial conditions used in the domain of attraction pro- 
cedure are found as follows. The variables x, through x , are set as 

1 n-1 

random variates uniformly distributed over their respective ranges. The 
final state variable is found by subtracting the sum of the n-1 vari- 
ables from the total mass. If x is within its limits, the n state vari- 

n 

ables are accepted as members of the set. If the last variable is out of 
its range, the entire collection of n variables is discarded. This pro- 
cess is repeated until the required number of initial conditions are 
found. 


The mass constraint forces the admissible initial conditions into 
nonuniform distributions. Figure 32 shows histograms of the marginal 
distribution of the edible and waste storage initial conditions for the 
three-chamber model. The six-chamber model is shown in Figure 33. The 
other state variables have similar distributions. Few large values are 
present because once a large value is selected there is a high probabil- 
ity that the total-mass constraint will be exceeded. 

Such nonuniform distributions must considered in determining the 
domain of attraction's volume and minimum radius. Clearly, the fraction 
of initial conditions in the domain is not proportional to the domain's 
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Figure 32: Sample Distribution Histogram; Three-Chamber Model 
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volume. However, if the state space is divided into moderately sized, 
identical cubes (actually they are 9-, 16-, or 30-dimensional hypercubes 
in these models), the initial conditions can be assumed to be uniformly 
distributed in each. The number of points in each cube gives the local 
density. Therefore, the volume of the domain of attraction is: 


V 


K 

2 

i=1 



(5.2) 


where ( N p ) i is the number of pass points and d^^ is the local density of 
points in cube i. The number of cubes in the state space is K. 


Estimating the domain’s boundary requires knowledge about the local 
point spacing. Its expected value is determined by taking the n root 
of local volume per point for an n-dimensional state space. The minimum 
radius is the distance to the nearest fail point minus the local point 
spacing. 


The state space, while bounded, is very large for these two models. 

The three-chamber model has a 9-dimensional state space with a volume of 
1 4 

8.9x10 . A 10,000 point sample of initial conditions would have an 

average spacing of 16.4. The six-chamber model has a 16-dimensional 
20 

volume of 3.0x10 ; its average point spacing for 10,000 points is 10.7. 

2 6 

The 30-dimension twelve-chamber model has a volume of 8.0x10 with an 
average spacing of 5.8 for a 10,000 point sample. As shown in the dis- 
tributions of Figures 32 and 33, the local density and point spacing can 
be expected to vary considerably throughout the state space. 

15.2 System Behavior and Storage Capacity 


The CELSS model behavior is strongly dependent on the selection of 
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storage tank capacities. Here, variations of the waste tank capacity are 
considered. Simulations are made to delineate the domain of attraction 
using a random selection of 10,000 initial conditions. These initial 
conditions are viewed as the point the system arrived at after a random 
component failure. It is assumed that the component has been fixed so 
the subsequent model evolution can be followed with the discrete-time 
simulation. Proportional control with gain k = 1 is used in the examples 
that follow (see Figure 29). System mass is set so there is enough for 
10-day waste and food buffers. 

When the three-chamber model is used with waste storage capacity 
larger than 82.4, all 10,000 initial conditions are attracted to the 
equilibrium point that generates food at the rate the crew consumes it 
(see Figure 25). This steady state is therefore globally attracting and 
stable. It should be noted that this attractor is an equilibrium point 
only in the discrete-time simulation used here. In a continuous-time 
simulation it would be a stable limit cycle. Such a steady state is 
called the pass steady state or the pass equilibrium point because a 
CELSS operating at this point meets its crew's food and environmental 
needs. 


When the waste storage tank is smaller than 70, all 10,000 initial 
conditions are attracted to equilibrium points that have edible harvests 
of zero. Each value for the waste capacity has a corresponding, unique 
equilibrium that is globally attracting. All of these points are called 
fail points. 

In the intermediate range of storage sizes, from 70 to 82.4, both a 
pass and a fail equilibrium point coexist. The fraction of the state 
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space attracted by each is shown in Figure 34. These volume fractions 
have errors of +0.005. As the tank size increases, the region attracted 
by the fail equilibrium shrinks while that of the pass increases. For 
waste capacities below 70 the pass point is unstable and, for storages 
larger than 82.4, the fail point is unstable. Thus, there is a bifurca- 
tion and reverse bifurcation in the stable points as the waste capacity 
is increased. There are no other attractors in the three-chamber system. 

One possible design goal is to maximize the domain of attraction 
about the pass equilibrium. If this domain covers the entire state space 
then the system is able to recover from all perturbations. However, it 
is desirable to have a system with small storage tanks and mass. Figure 
34 indicates that all initial conditions are attracted when the waste 
capacity is 82.4. A typical, continuous-time trajectory for this system 
is shown in Figure 35. here it can be seen there are no edible harvests 
for 100 days. Further, it takes the system 240 days to return to the 
steady state. This trajectory is stable, but it is not survivable. 

A performance criteria can be established such that all pass ini- 
tial conditions must settle to the pass equilibrium point in t £ 300 
days. Figure 36 compares the volume of this performance domain with the 
domain of attraction. To meet the performance limit over the entire 
state space, a waste storage capacity of 83 is needed. The minimum 
radius as a function of waste storage size for the two domains is ident- 
ical (see Figure 37). Thus, the extra volume in the stable region 
appears to be far from the equilibrium. Other performance criteria could 
be used such as a limited time without food, etc. 


An alternative approach to reducing storage capacities is to reduce 
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Figure 35: Long Transient Settling Time; 

Three-Chamber Model (Continuous-Time) 
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Figure 36: Pass Domain of Attraction and Performance; Three-Chamber Model 
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total system mass. For instance, if a 5-day buffer is used in the food 
and waste tanks, the total three-chamber system mass is reduced to 116. 
The domain sizes for the pass and fail equilibria are shown in Figure 
38. As with the 10-day buffer example (Figure 34), there is a bifurca- 
tion pattern between the pass and fail attractors. However, the pass 
equilibrium is globally attracting for waste capacities above 73. 
Therefore, in addition to reducing system mass by 20, the waste storage 
can be made 9.6 units smaller without changing the system* s stability. 
Smaller buffers, however, cause there to be longer transients when food 
is not available. 

The six-chamber model has all of the same dependencies on storage 
capacities and system mass as the three-chamber model. However, it has a 
more complicated bifurcation pattern as well as attractors that are not 
equilibrium points. 

Figure 39 shows the percent of the state space attracted to the 
pass and fail equilibrium points. The pass steady state attracts the 
entire space for waste capacities larger than 54. The fail equilibrium 
is globally attracting for storages smaller than 43.2. While the waste 
capacity needed for the pass domain to cover the entire space is smaller 
than that needed in the three-chamber model, many of the initial condi- 
tions still take a long time to settle. This is shown in Figure 40 where 
the performance domain has a settling time limit of 300 days. 

In the intermediate range of waste storage capacity, the equili- 
brium points go through a series of bifurcations that create higher- 
dimensional attractors. The motion on these attractors is bounded and 
may be quite complex. For example. Figure 41 shows a trajectory that, 
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Figure 38: Domains of Attraction with Total Mass=116; Three-Chamber Model 





Figure 39: Domains of Attraction; Six-Chamber Model 
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Figure 41: Limit Cycle Trajectory; 

Six-Chamber Model (Discrete-Time) 
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after an initial transient, becomes periodic. This figure shows the 
output of a discrete-time simulation where the output samples have been 
connected to form the trajectory. Thus, the edible harvest plot shows 
harvest levels oscillating in a repeating pattern. This limit cycle 
behavior is globally attracting for a waste capacity of 46. 

When the waste storage capacity is 53.4, two-thirds of the initial 
conditions are attracted to the discrete-time trajectory shown in Figure 
42. This motion is called chaotic. It is characterized by trajectories 
that never repeat, a continuous spectrum, and a Lyapunov exponent 
greater than zero [10]. The first Lyapunov exponent was found to be 0.3. 
The chaotic motion takes place on what is called a "strange attractor" 
[14]. 

While it is desirable to use the smallest waste capacity possible, 
it is not Acceptable to go into this region of higher-dimension attrac- 
tors. The limit cycle and chaotic motion shown in Figures 41 and 42 are 
not survivable because there are long periods with no edible harvest and 
no food in the storage. No stable attractors capable of providing a con- 
tinuous food supply were found in the six-chamber model, other than the 
pass equilibrium point. 

In a twelve-chamber model, limit cycle behaviors can be found that 
meet the crew's needs. This system has a pass equilibrium point that 
provides harvests of 5 units every 5 days. The total mass with 10-day 
food and waste buffers is 113.5. There are 30 state variables in this 
model. Figure 43 shows the fraction of the state space attracted to the 
equilibrium points and higher-dimension attractors. A waste capacity of 
48 guarantees global stability but, as was true in the other models, the 
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Figure 42: Chaotic Trajectory; 

Six-Chamber Model (Discrete-Time) 
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settling time may be long for some of the initial conditions. 

A waste capacity of 37.3 gives a globally attracting cycle (see 
Figure 44). There is always a 1 unit per day flow of food to the crew 
making this cycle survivable. Further, it attracts the entire state 
space with a waste capacity smaller than that which would support a pass 
equilibrium point. Ordinarily, cycles are engineered out of systems 
because they wear out components. However, biological systems are not 
hurt by constant oscillation; in fact this behavior may be more natural. 
Another benefit is that repairs and maintenance can be performed without 
disturbing the system every 110 days as the low point in the flow goes 
through the component. 

This cycle is aperiodic. It has "envelope" and "carrier" frequen- 
cies that do not have a common divisor. The attractor creating this type 
of trajectory can be pictured as motion on a torus [10], However, for 
the purposes of insuring that the crew's needs are met, the aperiodic 
component of this motion is not significant. 

The twelve-chamber model also has chaotic motion. Figure 45 shows 
the motion (discrete-time simulation) on a strange attractor that is 
globally attracting at a waste capacity of 37.6. While there are long 
periods where the crew's food needs are met, such occurrences are diffi- 
cult to predict. It should be noted, however, that there may be strange 
attractors that can insure the needed food flow. These would be as 
acceptable, from a survival point of view, as the cycle shown in Figure 
44. To demonstrate the complexity inherent in a chaotic trajectory. Fig- 
ure 46 shows the edible harvests for 13.7 years. 
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Figure 44: Oscillating Trajectory; 

Twelve- Chamber Model (Discrete-Time) 
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Figure 45: Chaotic Trajectory; 

Twelve-Chamber Model (Discrete-Time) 
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In summary, changes in waste storage capacity and the number of 
plant chambers significantly alters the system's behavior. All three 
models show a pattern of bifurcation and reverse bifurcation as the 
waste capacity is increased. The three-chamber model has only two 
equilibrium points while the six- and twelve-chamber models have two 
equilibrium points and a collection of higher-dimension attractors. As 
the waste capacity is reduced, the settling time of the pass initial 
conditions increases. In addition to permitting the use of smaller 
storage tanks, the six- and twelve-chamber models have oscillating, 
stable, globally attracting trajectories. These occur at waste capaci- 
ties smaller than those that permit surviving equilibrium point 
behavior. 

5_.Z Information Flow in a CELSS 

Since the CELSS system is a closed cycle, it is not apparent what 
piece of information will be most useful to a controller. To investi- 
gate the effectiveness of different information use, the three-chamber 
model was chosen with a waste capacity of 82. The design goal is to max- 
imize the volume and the minimum radius of the domain of attraction of 
the pass equilibrium point. The domain is delineated with 10,000 initial 
conditions. 

In the first example, the waste flow to the processor is set at the 
steady state value of 3 units/day. This no-control situation is the 
equivalent of an open-loop investigation in traditional process control. 
The domain of attraction for this case is 0.056 + 0.005 of the state 
space volume. The minimum radius is 1.5 + 1.1. 
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A proportional control is added to the system. Its general form is: 


u 


u 


steady 


t1 + k<i - y/y set n 


(5.3) 


where u is the waste flow to the processor, u s t eac |y is the steady state 

processor flow, and k is the proportional gain. The observed variable 

is y and the desired value for y is y . The quality of the control 

sez 

depends on the selection of the state variable y that is used in the 
controller feedback. 


Figure 47 shows some possible paths for feedback information flow 
in a CELSS. Three feedback configurations are shown, for which domain 
of attraction simulations are repeated with 100 proportional gains ran- 
domly selected from the range: 

0 < k < 5 (5.4) 

The gain that provides the largest minimum radius is the controller for 
that feedback system that can recover from the largest set of random 
state perturbations. 

When the Waste flow is set by comparing the level in the waste 
storage with the steady state, the largest domain of attraction volume 
is 0.099 + 0.005 of the state space and its minimum radius is 2.1 + 1.1. 
Its gain is k = 0.51. It should be noted that volume fractions from 
0.050 to 0.099 can be obtained by using gains in the range from 0 to 
1.5. A gain of zero is the no-control situation. Using the waste storage 
level as the control input has a marginal effect on the system’s ability 
to recover from random failures. 


The presence of inert matter in the plant chamber reduces edible 
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growth (see Figure 24). A controller can compensate by increasing the 
nutrient flow to the plants, which increases the edible fraction. 
Searching the range of gains with this feedback design gives a maximum 
domain of attraction volume of 0.180 + 0.005 of the state space. The 
minimum radius is 5.2 + 1.1. This control has a larger domain than the 
waste- feedback or no-control configurations. Although the gain for this 
domain is k = 1.4, there is a wide range of gains that give comparable 
results. Volumes from 0.100 to 0.180 are achieved with any gain larger 
than k =0.5. As the gain gets larger, the controller reaches its max- 
imum flow for very small inert masses. Such saturation indicates that 
most results above a certain gain will be identical. 

The harvest can also be used as input to the controller. Although 
this configuration is like the one used in the waste capacity investiga- 
tion of the previous section, here the capacity is fixed and the control 
gain is searched to find the system that recovers from the largest set 
of random failures. The largest domain of attraction was found for the 
gain k = 1.2. The domain’s volume is 0.585 +^0.005 of the state space 
and its minimum radius is 18.3 +1.1. The domain is rather sensitive to 
gain. For gains from 1.0 to 1.25 the domain volume fraction is between 
0.500 and 0.585. However, there is a sharp drop in domain volume for 
gains outside this range. 

Global system behavior was slightly improved over the no-control 
situation by using the waste storage level or inert mass as the feedback 
information. However, a large increase in the domain of attraction was 
obtained when harvest information was used. This can be understood by 
redrawing the CELSS in a more conventional process control configura- 
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tion. Figure 48 shows the system drawn with forward and backward paths. 
Control is exerted on the system by adjusting the waste processor flow. 
Using the waste storage level as the control input gives very little 
improvement in the size of the domain of attraction. Since this domain 
represents stability, it is expected that feedforward control would have 
no influence on its size. Hence, making use of the waste storage infor- 
mation mimics the traditional behavior of a feedforward path. The inert 
mass in the plant chambers provides feedback on the inner mass loop 
which gives more system stability. Feedback on the outside mass loop, 
using harvests as the control input, does even more to improve the sta- 
bility domain. 

This analogy with process control helps to explain the effective- 
ness of various feedback configurations. It is, however, only a single 
input, single output controller investigation. More complex controllers, 
such as those that use state feedback, would probably be able to achieve 
larger domains of attraction and performance than these output feedback 
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6. SUMMARY 

In this dissertation, a technique for analyzing and designing con- 
trollers for nonlinear systems based on measures of the domain of 
attraction has been presented. These methods are particularly suited to 
investigating the dynamic consequences of changes in the waste storage 
capacity, the system mass, and how information is used for control in 
CELSS models. The models' high-dimensionality and nonlinear state equa- 
tions make them difficult to analyze by any other technique. 

The domain is the region of initial conditions that asymptotically 
approach an attractor. The attractor may be an equilibrium point, a 
limit cycle, or other higher-dimension attractors. A refinement of this 
region is the domain of performance which is the region of initial con- 
ditions that meet a performance criteria. 

In nonlinear systems, local stability does not insure stability 
over a larger region. The domain of attraction marks out this stability 
region. In this way, it is a measure of nonlinear system's ability to 
recover from random state perturbations. In linear systems, local sta- 
bility guarantees global stability so this concept is not enlightening. 
However, the use of a domain of performance as a global measure is use- 
ful in both linear and nonlinear systems. 

When considering random perturbations, the minimum radius of the 
domain is a measure of the magnitude of perturbations for which recovery 
is guaranteed. An advantage of this measure is that it is a vector 
length and, therefore, is applicable to systems with any dimension. The 
representation of global system performance in a single scalar permits 
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easy comparisons of system and controller designs. 

The domain of attraction or performance is delineated by randomly 
selecting initial conditions from the region of state space being inves- 
tigated. If these are from uniform distributions, the number of points 
in the domain is proportional to its volume. The minimum radius of the 
domain’s boundary is found by subtracting the typical point spacing dis- 
tance from the nearest initial condition not in the domain. In con- 
strained systems, the sample distribution cannot be uniform. The minimum 
radius is now found by subtracting the local point spacing distance from 
the nearest point not in the domain. By numerically inverting the 
nonuniformity of the sample points, the domain volume can be obtained. 
Error determinations for the volume and minimum radius are found by 
repeating the domain simulations with independent samples. 

The use of the domains of attraction and performance were demon- 
strated in controller design for an inverted pendulum. This global 
design technique was contrasted with more traditional, local design 
methods. Since the problem is nonlinear, the traditional designs only 
insured some sort of local behavior. There is no map between local 
design decisions and the nonlinear system's global behavior. Thus, it is 
not apparent what parameter changes need to be made to improve the glo- 
bal performance. 

When the domain of attraction or performance is used in the con- 
troller design, the global behavior of the system is used as the design 
criteria. Consider the case where a controller form is selected and the 
gains need to be set. A controller can be designed by randomly selecting 
gains, simulating the resulting domains, and using the minimum radius as 
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the selection parameter. For the inverted pendulum, this technique was 
used to find a linear state-feedback controller that was able to recover 
from a larger class of state perturbations than those found by local 
techniques. A nonlinear controller was also designed which had larger 
domains of attraction and performance than the linear one. 

An advantage to this global design technique is that the design 
process is automated once the system model is formulated, the region of 
state space to be investigated is established, and the form of the con- 
troller is set. The domain simulations are made from randomly selected 
initial conditions over the state space. Each domain simulation results 
in a minimum radius representing the system's ability to recover from 
random state perturbations. The complexity of the controller has no 
effect on this process, other than to increase the number of parameters 
that need to be set for each domain simulation. The system dimension is 
not a problem because the selection is made using a scaler. 

The domain of attraction was used in the CELSS examples not only to 
find good control gains, but also to investigate system behavior as a 
function of system parameters. It was found that the three-chamber model 
has only two equilibrium points in the discrete-time simulation. One 
point provides edible harvests that meet the crew's food requirements. 
The other has edible harvests that are always zero. The former, surviv- 
ing equilibrium is called the pass equilibrium point and the latter one 
is called the fail equilibrium. For small waste storage capacities the 
fail point is the only attractor. As the waste storage increases this 
equilibrium bifurcates into pass and fail stable points. Then there is a 
reverse bifurcation leaving only the pass equilibrium point. 
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The six- and twelve-chamber models showed a similar bifurcation, 
then reverse bifurcation pattern. However, in the intermediate range of 
tank capacity, there were also stable, higher-dimension attractors 
present. These have bounded trajectories that are limit cycles, 
aperiodic oscillations, and chaotic motion. While there were no surviv- 
able oscillating trajectories found for the six-chamber model, the 
twelve-chamber system had survivable behaviors that are cyclical. 

An investigation of the effectiveness of using different informa- 
tion as controller input showed that the CELSS can be viewed in a pro- 
cess control configuration. Here, controls that look like traditional 
feedforward control have no effect on the system stability. Feedback 
control, however, improves system stability. In this representation the 
multi- loop form of the CELSS becomes apparent. 

To achieve the minimum mass and storage capacity in a CELSS, sys- 
tems with more frequent harvests are preferred. Less mass is tied up in 
the plant chambers and the harvests are small, requiring less storage 
capacity. However, these systems have a large number of state variables. 
Real-time controllers may be burdened by this high-dimensionality. 

Smaller waste capacities may be used to force the system into a 
cyclical mode. While this is not usually acceptable in mechanical sys- 
tems, it appears to have benefit in a CELSS. The waste capacity is 
reduced while not requiring an increase in the food storage capacity. 
Also, by waiting for the low point in the system flow, repairs and 
maintenance could be performed without disturbing the system. 


The abstract CELSS models presented in this dissertation only give 
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a hint of the complexity of dynamic behavior that will be observed in a 
real system. Modeling the atmosphere as a separate loop would probably 
create new system behaviors due to the interaction of the gaseous and 
solid/liquid loop. State feedback controllers may be successfully 
designed to globally stablize the system about a variety of operating 
behaviors. 
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